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Abstract 

The sedimentation of a heavy Stokes particle in a laminar plane or axisymmetric flow is in- 
vestigated by means of asymptotic methods. We focus on the occurrence of Stommel's retention 
zones, and on the splitting of their separatrices. The goal of this paper is to analyze under 
which conditions these retention zones can form, and under which conditions they can break 
5_j \ and induce chaotic particle settling. The terminal velocity of the particle in still fluid is of the 

order of the typical velocity of the flow, and the particle response time is much smaller than the 
^5 ■ typical flow time-scale. It is observed that if the flow is steady and has an upward streamline 

£S) , where the vertical velocity has a strict local maximum, then inertialess particle trajectories can 

take locally the form of elliptic Stommel cells, provided the particle terminal velocity is close 
,— ,■ enough to the local peak flow velocity. These structures only depend on the local flow topology 

and do not require the flow to have closed streamlines or stagnation points. If, in addition, the 
| flow is submitted to a weak time-periodic perturbation, classical expansions enable one to write 

the particle dynamics as a hamiltonian system with one degree of freedom, plus a perturbation 
containing both the dissipative terms of the particle motion equation and the flow unsteadiness. 
Melnikov's method therefore provides accurate criteria to predict the splitting of the separatrices 
of the elliptic cell mentioned above, leading to chaotic particle trapping and chaotic settling. 
The effect of particle inertia and flow unsteadiness on the occurrence of simple zeros in Mel- 
nikov's function is discussed. Particle motion in a plane cellular flow and in a vertical pipe is 
O I then investigated to illustrate these results. 

— , Key-words : Two-phase flows, Stokes particles, sedimentation, chaotic motion, Melnikov's 

method. 
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1 Introduction 

The motion of tiny heavy particles in either laminar or turbulent flows is of major importance for 
the understanding of many industrial or natural flows. Even when particles do not significantly 
modify the surrounding flow, this topic shows many unsolved questions. When the flow is 
turbulent, modeling the dispersion and sedimentation of heavy particles is a difficult task which 
has retained much attention since the pioneering works by Csanady pQ, Snyder & Lumley [2] 
or Maxey [3], to name but a few. In laminar flows, particle motion is also investigated for at 
least two reasons. First, a better understanding of particle dispersion and sedimentation in fully 
developped turbulence can be achieved if one understands the detailed interactions between 
particles and elementary flow structures like vortices or shear layers. Secondly, even in laminar 
flows, the behaviour of small heavy particles can be rather unexpected, and chaotic settling 
or permanent suspension ("trapping") can occur. In the limit of vanishing particle Reynolds 
number, these complex trajectories are not due to wake effects such as vortex shedding, since 
the disturbance flow set up by the inclusion obeys the linear Stokes equation. They are due 
to the gradients of the unperturbed fluid flow, which introduce severe non-linearities into the 
particle motion equations. 

Experiments and calculations conducted since 1900 revealed that passive heavy Stokes par- 
ticles dropped in a cellular flow can be trapped in the vicinity of vertical upward streamlines 
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and remain suspended for a long time, under the combined effect of gravity and of the hydrody- 
namic force (Benard [3], Stommel [5], Simon & Pomeau [6J, Cerisier aZ. [7J, Tuval ei ctZ. [5]). 
Such nearly closed trajectories are often called "Stommel retention zones", even though they 
are rather loosely defined. They are the central point of the present work. Stommel investigated 
these structures to understand plankton settling in a lake, where wind- induced cellular motion 
is likely to induce strong inhomogeneities in plankton concentration. In the last decades, various 
numerical analyses revealed that aerosol sedimentation can display complex features, whether 
chaotic or not, in elementary laminar flows such as plane cellular flow (Maxey & Corrsin [9], 
Fung |10j . Rubin et al. |llj). or ABC flow (Mac Laughlin [12J). Some of these analyses revealed 
that the occurrence of Stommel cells in various parts of the flow played a key role in complex 
particle motions. For example, some particles released in a plane cellular flow, submitted to a 
weak time-periodic perturbation, were observed to wander around some kind of Stommel cells, 
and could also be trapped temporarily within these cells, in a rather chaotic manner [10]. The 
occurrence of Stommel cells, as well as their role in chaotic sedimentation, is therefore a topic 
of interest to understand and quantify complex particle motion in simple flows. 

The goal of the present paper is to derive analytical criteria giving the particle parameters for 
which chaotic sedimentation and/or trapping occurs, for heavy isolated Stokes particles dropped 
in the vicinity of a vertical streamline where Stommel retention zones are likely to form. Two- 
dimensional and axisymmetric divergence-free flows will be considered. The paper is organized 
as follows. First, particle dynamics in the vicinity of an upward streamline of any steady 2D flow 
is investigated (section [2]) . Conditions leading to the formation of retention zones are derived. 
Then the splitting of the separatrices of these retention zones is analyzed by means of Melnikov's 
method (section [3|) . This analysis is applied to chaotic sedimentation in a plane cellular flow in 
section 01 Finally, these results are generalized to axisymmetric flows ( section E]) . In the next 
few paragraphs we present the basic equations of particle dynamics used in this paper. 

If X p (t) is the particle position at time t, the simplest motion equation for heavy isolated 
Stokes particles in a fluid of infinite extent is 

m p X p = m p g + 6n(j,f &(Vf(X p , t) - X p ^j , 

where m p and a denote the particle mass and radius respectively, ///is the fluid dynamical 
viscosity, and Vf is the fluid velocity field in the absence of particle. This equation is valid 
provided the particle Reynolds, Stokes and Taylor numbers are much smaller than unity : 

Re = J—E 11 <^ i an d St = < 1, 

v v 

and Ta = — 11 < 1, 

v 

so that the disturbance flow due to the inclusion is a creeping quasi-steady flow (u?i is the 
typical time-scale of this flow). In addition, the particle density p p must be much larger than 
the fluid density pf , so that buoyancy and pressure gradient force of the undisturbed flow can 
be neglected. Also, brownian motion is assumed to be negligible. The velocity field Vf is taken 
to be a steady flow submitted to a weak unsteady perturbation : 

V f (x,t) = V f 3 (x)+eV f 1 (x,t), e<l. 

Following Fung |10j we will consider the case where Vj is T-periodic with T = 2-k/uj. If Vo 
denotes a typical fluid velocity and Lq is the length scale of velocity gradients in the absence of 
inclusion the particle motion equation, non-dimensionalized by (Vq,Lq), is : 

tX p = V t + Vf{X p ) + e Vj{X p , t) - X p (1) 
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without renaming t, X p and Vf. The parameter r is the non-dimensional response time, which 
is assumed to be much smaller than unity : 

-> -> 2/9 a?g 1 

and Vt is the non-dimensional terminal velocity in a fluid at rest ( Vt = — — ), which is 

9p f v V 

assumed to satisfy : 

V T = \V T \ = 0(1). 

The Proude number t/Vt = Vy/gLo is therefore much smaller than unity. Also, the terminal 
particle Reynolds number is of order (a/Lo) VoLq/v, and is therefore very small provided the 
flow Reynolds number VqLq/^ is not too large. The particle motion equation contains two 
independent small parameters, namely e and r. In the following we set r = ke, where k 
denotes any parameter held fixed as e — > 0. Then classical asymptotic expansions [3] of the form 

X p = V T + Vf + eVp 1 + 0(e 2 ) lead to : 

X p = V T + Vf(X p ) 

+ e(v f 1 (X p , t) - kWf.{Vf + V T )) + 0{e 2 ). (2) 

The case k = corresponds to heavy particles without inertia [9][3][6j : their velocity is equal 
to the fluid local velocity plus Vt- In this case the particle dynamics is non-dissipative since 
div Vf = 0. In the following we will consider the general case k > 0. 



2 Formation of retention zones in plane flow 

As noticed by many authors [5 J [9 J [6 J [lOj, when the flow is two-dimensional the dynamical system 
([2]) is, to leading order, hamiltonian : 



*p=( • )=V T + Vf(X p )=( % ) (3) 



8H 
dy 

m 

dx 



where H(x,y) = ip°(x,y) + xVp, and Tp° is the streamfunction associated to V9, and x is 
the horizontal coordinate and y is the upward vertical coordinate. This enables to notice the 
following property. Suppose ip°(x,y) is of class C°° and has a vertical upward streamline S (Fig. 
[T](a)) and a point C € S where the vertical velocity v = —dip°/dx satisfies : 

dv , s dv , , 

d^v d^v / d^v \ 2 

dx 2 ^ ^ dy 2 ^ ^ \dxdy^ ) ^ 

d 2 v 

jpOQ) < o. (6) 

These are the classical conditions for v(x, y) to have a strict local maximum at C. We will denote 
x c and y c the coordinates of C, and x c = (x c ,y c ) its position vector. Also, we write V c = v(C) 
(the local peak flow velocity) and assume V c > 0. Then for Vr < V c and Vt sufficiently close to 
V c , the particle trajectory at order e° (Eq. (J3j) ) takes the form of an elliptic dipole of centre C 
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(Fig. QJb)). If in addition vorticity is zero along S , the axes of the ellipse are (C, x) and (C,y) 
and the semi-axes are 



I tyVc ~ Vt) , , = 2(V C - V T ) ( 

^(C) V ^V(C) 1 ' 



,xxx\^J V T ,xyy\ 

(where the coma indicates spatial differentiation). 



(a) 



y(t) 



(b) 



H(x,y) = constant 



2b 



2a E, 




x(t) 



Figure 1: Sketch of the local streamlines (a) of the unperturbed flow V9, and of the leading 
order phase portrait of particle dynamics (b) in this flow, under the hypotheses listed in the 
text. 



This property can be shown by expanding H (x, y) in the vicinity of C and making use of the fact 
that many derivatives of ip° vanish at C. First, u(x c , y) = for all y £ S , so d n u/dy n (x c , y) = 
for all integer n, so d n t^° /dy n {G) is zero too for all n. Also, velocity is divergence- free, so 
ifj° xy (C) = u >x (C) = u x (C) + u, y (C) = 0. Also, ip° xx (C) = -v, x (C) = 0. Therefore, by expanding 
H(x, y) we are led to : 

ii ... , - ^ 



H(x,y) = ^{G) + x c V T + {x 
Vt-Vc + 



i .xxx / \2 , t xxii i \i \ 

(x - x c ){y - y c ) 



6 

+ 



I X Xr 



+ 



2 

+ K 



(8) 



where all the derivatives of have been taken at C and 1Z = 0(\x — x c \ A ). Therefore, the 
isolines H(x, y) = constant take locally the form of elementary curves. In particular, if 1Z is null 
or negligible, the line H(x,y) = H(C), for x / x c , is the quadratic curve : 



q(x - x c f + 2/3(x - x c )(y - y c ) + j(y - y c f + 5 = 



where 



a 



^xxxi^) n Ip^xxyiC) __ Ip^xyyiC) 



P 
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6 4 7 ' 2 

and 5 = Vt — V c - A sufficient condition for this curve to be an ellipse is : 5 ^ and 8 (a + 7) < 
and cry > j3 2 . Assumption (JUJ) implies a > 0, and by (0) we obtain 7 > 0. Also, Vt < V c means 
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S < 0. The two first ellipticity conditions are therefore fulfilled. The last one is a consequence 
of (J3J) . Indeed, Eq. (0) can be re- written as (6a) (27) > (4/3) 2 . Therefore we have 07 > §/3 2 , 
hence cr/ > (3 2 . 

If in addition vorticity is zero along S , we have tp° xx + ip° yy = so ip° xx = at any point in S . 
So tp° xxy is also zero along S . In particular fi = 0, and the ellipse has horizontal and vertical 
axes. In this case the semi-axes of the ellipse are a = \J ' —5/ot and b = \J— J/7, and this leads 
to (0). 

The remainder 7£ is negligible in the limit where \x — x c \ <C 1. On the ellipse this condition 
is fulfilled if max(a, b) -C 1, that is | Vt — V c |<C 1, since both a and 7 are 0(1). Therefore, 
the elliptical structure is asymptotically valid if Vt is sufficiently close to V c . In the following 
we will assume that 1Z is indeed zero or negligible, for the sake of simplicity. An analysis of the 
effect of 7Z is presented in Appendix O 

Under these conditions the particle dynamics in the vicinity of C has two stagnation points 
A = (x c , y c + b) and B = (x c , y c — b) located along the vertical direction, related by three 
separatrices (heteroclinic trajectories). These cells are an example of Stommel's retention zones, 
already observed with various shapes in many experimental or theoretical works cited in the 
introduction [lj [5] [6] [10] [7] [8]. To leading order e°, particles which are located within 
the elliptic cell are trapped and cannot exit. Also, particles released outside the cell cannot 
penetrate into the cell : these particles drop. Therefore, if one injects particles at random in the 
two-dimensional flow with a uniform distribution, the percentage of trapped particles should be 
approximated by the volume fraction of the elliptic cells, that is : 

= Ne rr a b = == = (Vc ~ V T ) (9) 

where _/V e is the number of elliptic cells per unit volume. (Since the problem is two-dimensional, 
volumes must be understood as surfaces multiplied by the unit length along the z direction.) 
This formula will be checked below (section d|) . 

If the 0(e) perturbation appearing in Eq. ([2]) is taken into account, the curved separatrix 
can be broken : particles could therefore go in and out in a quite unpredictable manner. The 
calculation of the parameters leading to separatrix splitting is therefore an interesting topic for 
the understanding of particle sedimentation, and is done in the next section. 

3 Heteroclinic bifurcation and chaotic particle motion 

For plane flows the particle motion equation ([2]) has the form of a hamiltonian system with 
"one-and-a-half" degree of freedom. The 0(e°) terms contain the fluid velocity as well as the 
sedimentation velocity, the combined effect of which can lead to the formation of the elliptic cells 
described above. The separatrices of these cells form heteroclinic cycles between the stagnation 
points A and B. The 0(e) terms manifest the effect of both the flow unsteadiness and particle 
inertia : clearly, these physical mechanisms might influence the particle dynamics in the vicinity 
of the elliptic cells, since phase portraits of the form of Fig. \H[b) can be structurally unstable. 

In order to investigate analytically the role of the 0(e) terms, we make use of the classical 
Melnikov method [13] to predict the appearance of intersection points between the stable and 
unstable manifolds associated to the hyperbolic points A' and B' of the Poincare section (x(to + 
nT),y(to + nT)), n G Z, of the perturbed system (Fig. [2]). (T is the period of the perturbation 
velocity field : Vj(x,t + T) = Vj(x) for all x and t.) Indeed, since the leading-order particle 
dynamics displays two hyperbolic saddle points A and B, the Poincare section of this system 
will have hyperbolic points too, at the same position A and B. If e is small enough the Poincare 
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B'(e) 



*(f„ +nTJ 



Figure 2: Sketch of the Poincare section of the particle dynamics at order e. In case (a) the 
manifolds do not intersect. Particles released from above the cell will not be trapped. Particles 
released within the cell will spiral out and exit the cell. In case (b) the manifolds intersect and 
chaotic sedimentation occurs. 



section of the perturbed system will therefore have two stagnation points of the same type 
(hyperbolic-saddle) A'(e) and B'(e), continuous functions of e, by virtue of the implicit functions 
theorem. Accordingly, there exists stable and unstable manifolds W s and W u attached to B'(e) 
and A'{e) respectively. If an unstable manifold intersects a stable one (e.g. point I of Fig. 12(b)), 
an infinity of such intersection points will exist, since both W s and W" are invariant subsets 
of the Poincare application. Because the O(s ) phase portrait displays a heteroclinic cycle, the 
particle dynamics in the vicinity of the broken separatrix corresponds to a horse-shoe map, by 
virtue of the Smale-Birkhoff theorem, characterized by chaotic trajectories and sensitivity to 
initial conditions. Particles are therefore likely to penetrate into the cell or exit from it, in a 
chaotic manner. 

The occurrence of intersection points can be predicted from the asymptotic calculation of the 
"algebraic distance" d(to) between W s and W u . It is defined as d(to) = SU.N, where N is the 
unit vector perpendicular to V} + Vr at a given (arbitrary) point P of the separatrix, and such 

that (V} + VT,N,e z ) is right-handed, and S (resp. U) are the intersection points between the 

line (P,N) and W s (resp. W"). (e z is the direct unit vector in the direction perpendicular to 
(x,y).) These points are shown, for the separatrix Si, in Fig. [2(a). At order e this distance is 
proportional to the Melnikov integral [13J : 

f + OO r 

M(t 



I +co r - -> i r - / - 

/ Vt + V} _ A V}(x (t),t + t 

J-oo 1 J ^o(*) L v 



kW}.[V} + V T 



X (t) 



dt (10) 



which is integrated along the separatrix of the unperturbed system, parametrized by Xq (t) with 
Xq(-oo) = A and Xq(+oo) = B and Xq(Q) = P. (Here, u A v is the triple product e z .(u A v).) 

In particular we have X = V}(X ) + V T and X = VV}.(V} + V T ). If particle inertia is taken 
into account (k ^ 0) the perturbation is dissipative, but this does not alter the validity of Eq. 
(|10p which only requires the unperturbed system to be non-dissipative. The Melnikov integral 
therefore reads : 

M(t )= X (t)AVUx (t),t + t )dt-km (11) 

J —oo 

where m = m(Vr) is independent of t$ : 

r+oo -, 

m{V T )= / X Q {t)AX {t)dt. 



This term is also independent of Vj and corresponds to the effect of particle inertia. Because 

Xq AXq is proportional to the curvature of the separatrix at Xo(t), the constant m manifests the 
contribution of centrifugal effects. One can easily check that m(Vr) < (hence —km(Vr) > 0) 
for the separatrix Si of Fig. QJb), and m(Vr) > (hence — km(Vr) < 0) for S3. This means 
that particle inertia tends to maintain SU in the direction of N on separatrix Si (that is SU 
points outward), and in the direction of — N on separatrix S3 (that is, SU points outward also 
since N points inward there, as (ty + Vr,N,e z ) is right-handed). In both cases this means 
that inertia tends to maintain W u at the periphery of the unperturbed cell and W s towards the 
interior of the cell. In particular, this enables to conclude that in the steady case (u) = 0) the 
phase portrait of the 0(e) particle dynamics takes the form of Fig. [2£ a). The Stommel cell is 
broken because of particle inertia only : particles released above the cell will never be trapped, 
and particles released inside the cell will always escape : these results agree with those of Rubin 
et al. [11], which have been obtained for aerosols in a cellular flow field. 

In the following, the Melnikov integral is calculated in the case where 1Z = in Eq. and 
for the elliptic separatrix Si of Fig. QJb). Also, the constant (3 is taken to be zero, for the sake 
of simplicity. The calculation of Xq{€) is shown in Appendix [A] and we are led to : 



We observe that m(Vr) is proportional to the ratio a/b (if Vt — V c is held fixed), so that flat 
horizontal ellipses (b <C a) will be more sensitive to particle inertia than flat vertical ellipses 
(b S> a), as expected for a centrifugal effect. 

The non-constant part of the Melnikov function requires to specify the kind of perturbation 
applied. In this paper, "homothetic" perturbations of the form (Fung |10j ) 

Vf(x,t) = Vf(x) sinwt, (13) 

will be considered. This corresponds to flows with variable intensity without any change in 
the shape of the streamlines. Clearly, there is no chaotic advection of pure tracers in this 
flow, but it will be shown below that chaotic particle sedimentation will occur. In the case 
of a homothetic perturbation the term V® A Vj vanishes, and calculations lead to M(to) = 
— A(ui, Vr)cosujtQ — km(Vr), where A(uj,Vt) is proportional to the Fourier transform of the 
horizontal component of Xo(t) and reads (see Appendix |A} : 

IM VS itVtuj 
A(w, V T ) = , ^ nuJ , (14) 

\J ^%yy txxx C ° Sa 2V2^° xyy (V c -V T ) 

where derivatives of have been taken at C. Clearly, if A(cj, Vt) < k \ m(Vr) \ the Melnikov 
integral does not have any zero, and this provides an analytical criterion to predict the splitting 
of the curved separatrix of the retention zone. It also enables to predict the lowest velocity Vt 
below which no chaotic trapping can occur (see appendix IB"]). Note again that these calculations 
are expected to be valid either if 1Z = or if Vt is sufficiently close to V c . These results will be 
used in section [5] and compared to numerical solutions of Eq. ([1]) . 



4 Application to plane flows 

In order to validate the above results we have computed numerical particle trajectories for 
various plane flows, and checked whether chaotic settling appears for the parameters predicted 
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by Melnikov's analysis. The results are particularly accurate when the basic flow is itself a third- 
order polynomial function of (x, y), e.g. for a cellular flow of the form ip°(x, y) = x(x 2 + y 2 — 1), 
since the Taylor expansion ([8]) is exact there. Results are also satisfactory for various non- 
polynomial flows, provided | Vr — V c \ is small enough. Consider for example the streamfunction 
ijp(x,y) = sinx siny, already investigated by several authors (Maxey & Corrsin [9], Fung [TO] . 
Rubin et al. [11]). One can easily check that there is an array of upward vertical streamlines in 
the unperturbed velocity field, for example : x = and —it < y < 0. The vorticity is zero on 
this streamline and the upward velocity has a maximum at C = (0, — 7r/2). Also, one can easily 
check that Eqs. are satisfied at C = (0, -vr/2), and that V c = 1. So if V T < 1, and 

Vt close to 1, the particle trajectories at order e° in the vicinity of C take the form of elliptic 
dipolar vortices with half-axes given by Eq. ([7]) : 

a = y/6(l - V T ) and b = ^2(1 - V T ). 

Figure [3|a) shows typical numerical trajectories in this flow when Vt = 0.8, in the steady 
case and with very small particle inertia (r = 0.002). The Stommel cell is clearly visible, in 
agreement with the numerical computations of Refs. [5] [9] [10]. The ellipse of semi-axes a and 
b is also plotted for comparison. In order to check the theoretical result Q, which is a direct 
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X 

Figure 3: Upper graph : typical particle trajectories when Vt = 0.8, to = 0, e = 0.001 and 
k = 2 (thick solid lines, numerical computation). Thin solid lines are Bow streamlines. The 
dashed line is the ellipse corresponding to the analytical result (0). Some particles have been 
released above the cell, others have been released inside the cell. Lower graph : trajectories of 
two particles released slightly above the cell when e = 0.03, r = 0.06 and oj = 0.5. 

consequence of the appearance of the elliptic cells, we have uniformly released 5000 particles 
within the domain [— n, n] 2 and solved numerically the full motion equation (pQ) in the steady case 
(oj = 0), and with a finite but small particle response time r. Then particles above the horizontal 
line y = —it at time t have been called "suspended" . The percentage of suspended particles for 
large t is expected to be close to the volume fraction of the elliptic cells $ v . Figure H] shows this 
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percentage versus t, together with & v (Eq. ([9]) with N e = 2/(2ir) 2 since there are two elliptic 
cells in each periodic box), for Vt = 0.8, and in the cases r = 0.001 and r = 0.05. Clearly, the 
percentage of suspended particles is close to for long times, but keeps on decaying because of 
particle inertia effects. The decay is almost undistinguishable for r = 0.001 but is clearly visible 
for r = 0.05. Indeed, since r is finite in the numerical solution, the Poincare section of particle 
dynamics is of the form of figure Eta), since the Melnikov function M(to) = —km is constant 
and non-zero in the steady case, so that the elliptic cells are not completely closed and tend to 
slowly empty under the effect of particle inertia alone (see also Rubin et aJ. [llj). 
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Figure 4: Decay of the percentage of suspended particles in a steady cellular Bow (numerical 
solution of Eq. {I])), with Vr = 0.8. The solid line is the analytical result (0) corresponding to 
t = 0. 

In the unsteady case (to > 0) the diagram (u),Vt) indicating the occurrence of separatrix 
splitting is obtained by plotting the curve A(lu,Vt) — \km(Vr)\ = where A(uj,Vt) is the 
amplitude of the Melnikov function (j!4j) and km is its constant part. This diagram is shown 
in Fig. [5] in the case (k = 2). If (oj,Vt) is located in the zone A(co,Vr) — \km(Vr)\ > then 
the separatrix is broken and particles can penetrate into the elliptic cell in an unpredictable 
manner : "chaotic trapping" occurs. The lowest terminal velocity below which no chaotic 
trapping is expected to occur is given by the analytical expression obtained in Appendix [Bj and 
is approximately Vt ~ 0.76. 

In order to compare these predictions of separatrix splitting with numerical solutions we 
have solved Eq. ([1]) with an implicit Runge-Kutta algorithm, for 500 particles released at t = 
above the elliptic cell. For t > these particles drop and go round the elliptic cell, but some 
of them are likely to penetrate inside in case of separatrix splitting, and this will significantly 
increase the length of their path. The computation of each inclusion is stopped if the particle 
reaches the "ground" y = —tt or if t = tf S> 1 (here, tf = 50). We then measure the relative 
normalized averaged particle path length : 



C{u>) - £(0) 

W) ■ 



(15) 



where C(oj) is the length of the particle path averaged over the particles. This quantity is 
good indicator of the occurrence of separatrix splitting here. Indeed, if no trapping occurs 
then 9(uj) ~ 0, because all particle paths are rather similar (they look like those of figure 
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Figure 5: Plot of the bifurcation diagram given by Melnikov's analysis for sedimentation in the 
plane cellular Bow. 



[3(a)). If trapping occurs, some particles will spin inside the cell for a while (like in figure 
[3(b)) and this will increase 9{ui). Fig. [6] shows 9{oj) versus to (black dots), together with 
A(uj, Vt) — \km(Vr)\ (solid line), the positive values of which imply separatrix splitting (here 
e = 0.03 and r = 0.06). Clearly, when A(lo, Vp) — \km(Vr)\ > some particles are trapped, 
and when A(lu, Vp) — \km(VT)\ < no trapping is observed. This confirms that heteroclinic 
bifurcations do play an important role in the particle sedimentation process considered here. 



Note that for this flow one can easily check that TZ = 0(\ x 



') instead of 0(| x 



x r 



and it is shown in Appendix ([C]) that this allows us to neglect TZ provided 
enough. 



Vp — V c I is small 



5 Generalization to axisymmetric flows 

The analysis presented above can be readily generalized to particle sedimentation in an ax- 
isymmetric flow displaying an upward streamline S along the symmetry axis. If z denotes the 
coordinate along S and r is the distance to S , the velocity field can be written : 



V?.e r 



Vf.e z 



13V>° 
r dz 



1 dtp' 



o 



r dr ' 



where is the streamfunction of the unperturbed velocity field, which is assumed to satisfy 



ip°(r,z) = r 2 f(r,z) 



(16) 



where all the derivatives of / remain finite as r — > 0. The particle velocity in cylindrical 

coordinates reads X p = (r,z), and from Eq. ([2]) we obtain the particle dynamics at order e° : 

1 di{j p 



r = — — — and z 

r oz 



r dr 



where 



1. 



ip p (r, z) = ip°(r, z) + -V T r 
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Figure 6: Plot of the numerical computation of the averaged particle path length in the unsteady 
plane cellular flow tp = sinxsiny (black dots), together with A(ui,Vt) — k\m(Vr)\ (solid line, 
analytical expressions <fl2j) and ifH]) ). The fluid velocity at C is V c = 1, and r = 0.06, e = 0.03, 
(k = 2). 
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is the "particle streamfunction" . Suppose there exists a point C (r = 0,z = z c ) in S such that 
V9.e z has a strict local maximum. Then, in the vicinity of C, a calculation similar to the one of 
section [2] leads to 

2 

M r > z ) - y ( 5 + ar2 + - z c) + 7(2 - ^ c ) 2 ) 

plus higher order terms, with a = / )7T (C), j3 = f trz (C), 7 = f, zz (C) and 5 = Vr — V c . Here 
V c = V?.e z (C) = — 2/(C). The conditions ensuring that the the gradient of V?.e z vanishes at C 
and that the hessian is definite and negative there, are /, r (C) = /,z(C) = and 

8a7>9/3 2 and a > 0. 

Like in section [2l one can easily check that these conditions, together with Vr < V c , imply the 
ellipticity conditions 5^0 and 5(a + 7) < and 07 > /3 2 . The set ip p = is therefore an 
ellipsoid, together with the axis S . If, in addition, the curl of Vj is equal to zero along S , then 
one obtains (3 = 0, and the ellipsoid has vertical and horizontal axes. 

It is convenient to write the particle dynamics, in the vicinity of C, in a hamiltonian form. 
This can be readily achieved (in the case (3 = for the sake of simplicity) by setting p(t) = r 2 (t), 
so that 

< *» 



z 



' dp 



with H(p,z) = p(s + ap + j(z — z c ) 2 ^j + 7Z, where 1Z contains higher-order-terms. One can 
easily check that the phase portrait H = constant, when 5 < and 1Z = 0, is a closed cell 
bounded by a parabolic separatrix Si and a straight separatrix £2. The Melnikov function of 
the separatrix Si, when the flow is submitted to a homothetic perturbation (Eq. (|13|) ). can 
be calculated by using the same methods as above. After some algebra we obtain M(to) = 
—A(oj, Vt) cos(wto) — km(Vr), with : 

ttVt ^ 

^ = /rr(C)/«(C)smh( = - ) (17) 

and 



™i V T) = -\^J^{V c -V T fl\ (18) 

(Here also the comma indicates spatial derivation.) Like for the plane case, the diagram (w, Vt) 
indicating the occurrence of separatrix splitting is obtained by plotting the isolines of the curve 
A(u},Vt) — |A;m(Vr)|- This can be done as soon as the basic flow V9 is given, i.e. /, r r(C) and 
f,zz(C) are determined. These results are applied to a pipe flow in the next section. 



6 Application to particle settling in a vertical straight con- 
stricted duct 

Constricted ducts are investigated in various areas of fluid mechanics. For example, the axisym- 
metric creeping flow through a cylindrical duct with varying cross-section provides an elemen- 
tary biomechanical model for stenosed arteria or airways. It is known that the fluid acceleration 
through the constricted section can play an important role in various transfer processes, like 
particle motion. Here we consider a vertical straight constricted duct (Fig. [7]), with a vertical 
prescribed flow rate Q, laden with heavy microparticles released above the throat of the pipe. 
Under the combined effect of gravity and of the vertical flow, particles are likely to drop through 
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Figure 7: Sketch of the vertical constricted pipe. 



the throat, or drop to the throat and get deviated towards the pipe wall, or get trapped into 
the throat if the flow is perturbed. The non-dimensional streamfunction is taken to be : 

where R(z) denotes the pipe radius : -R(O) = R\ < 1 at the throat and R(oo) = 1 far from 
the throat. It corresponds to a parabolic Poiseuille flow in each cross-section of the pipe, with 
a prescribed flow rate Q = tt/2 and a peak axial velocity (at the tip of the parabola) equal 
to 1 far from the throat. This solution is valid if both the pipe Reynolds number and 1 — R± 
are small. We have chosen R(z) = 1 - (1 - Ri) exp(-z 2 ), with Ri = 0.8. At C = (0,0) we 
have f )r = f jZ = , f )ZZ > 0, / jTT > and f jZr = 0. The conditions required for V® .e z to have 
a strict local maximum at C are therefore satisfied, and we conclude that an elliptic Stommel 
retention zone can form there for particles such that Vr < V c = 1/Ri — 1.56 and Vp close 
enough to V c . The complete bifurcation diagram of this flow is obtained by plotting the curve 
A(u, Vr) — \ km(Vr)\ = 0, and is shown in Fig. [8] in the case k = 2. 
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Figure 8: Bifurcation diagram given by Melnikov's analysis for sedimentation in the constricted 
pipe how. 

These analytical results are compared to numerical solutions of Eq. ([T]) in Fig. for r = 2e = 
0.04 and Vr = 1-3. Like for the plane case we observe that the absence of particle trapping closely 
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Figure 9: Plot of the numerical computation of the averaged particle path length in the unsteady 
constricted pipe flow (black dots), together with the analytical expression of A(u, Vr) — k \m(Vr)\ 
(solid line, Eqs. Q2P and lflg[) ). The fluid velocity at C is V c ~ 1.56, and V T = 1.3, e = 0.02, 
r = 0.04. The absence of particle trapping corresponds to the absence of zeros in the Melnikov 
function (i.e. A(uj,V t ) - k \m(V T )\ < 0). 




2.5 



corresponds to the absence of zeros in the Melnikov function (i.e. A(uj,Vt) — \km(Vr)\ < 0). 
This agreement is observed even though the remainder 1Z of the Taylor expansion of H (p, z) is 
not strictly zero here. (Note that, like for the cellular flow investigated above, the remainder 
TZ(p, z) of this model pipe flow is of order 5 instead of 4.) 

In order to visualize the topology of chaotic particle dynamics we have plotted a particle 
cloud released at t = at the center of the throat (Fig. [T0j) . for the same parameters and with 
uj = 0.5. Trajectories are obtained by solving numerically Eq. (TjQ). The stretching and folding 
of the cloud is clearly visible, as expected for chaotic motion. 



t= T 



t= 2.5 T 



t = 



Figure 10: Plot of a particle cloud released at the center of the throat, for the conditions of 
figure and uj = 0.5. 



7 Conclusion 

The conclusions of this paper are valid for passive isolated tiny heavy low-Re particles with a 
response time much smaller than the flow time scales, and a terminal velocity of the order of the 
flow velocity. For such inclusions, sedimentation effects appear in the leading-order dynamics, 
whereas particle inertia can be treated as a perturbation, together with the flow unsteadiness. 
This enables to write the particle dynamics as a non-dissipative system submitted to a dissipative 
perturbation. 

Elliptic Stommel cells in steady plane or axisymmetric flows have been shown to occur as soon 
as the basic flow has an upward streamline displaying a strict local maximum with peak velocity 
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V c , and for heavy Stokes particles such that Vr < V c and Vr close enough to V c , in the limit of 
vanishing particle inertia. These results are not really new, as many authors already observed 
particle trapping in the vicinity of upward streamlines [5] [6] [7] [8]. However, the calculations 
of the present paper confirm that elliptic Stommel cells can form in a wide variety of flows, 
and do not require the flow to have closed streamlines or stagnation points. In particular, they 
can appear also if the basic flow is an open flow, like the vertical pipe flow with varying radius 
investigated above. 

When the flow is submitted to a weak time-periodic perturbation, the separatrices of these 
cells can be broken (heteroclinic bifurcation) and chaotic trapping, as well as chaotic trajectories, 
can occur. The parameters leading to separatrix splitting have been calculated in the case where 
only the quadratic non-linearities of the fluid velocity field can be taken into account (i.e. 7Z 
is zero or negligible). We have shown that particle inertia prevents the appearance of chaotic 
dynamics in the present case, by inducing an additive constant —km into the Melnikov function 
M(to). Therefore, particle inertia tends to maintain the stable and unstable manifolds W s and 
W u far from each other. This is clearly a centrifugal effect, since the constant part of the 
Melnikov function is mainly an integral of the curvature of the separatrix. In contrast, the 
flow unsteadiness tends to make the manifolds intersect, and the most dangerous perturbation 
frequency is the peak frequency of the Fourier transform of a combination of the coordinates of 
the point Xo(t) running on the separatrix. This combination depends on the detailed shape of 
the perturbation. The analytical expression of M(to) provides efficient criteria to predict the 
occurrence of chaos. These criteria have been tested by means of numerical computations, and 
the agreement is satisfactory : when ip° is a cubic function of position (1Z = 0) chaotic trapping 
is observed only if the Melnikov analysis predicts separatrix splitting. The agreement is also 
satisfactory if 1Z is non-zero but decays faster than | x — x c | 5 provided | Vr — V c | is small 
enough. 

The Melnikov function Mifo) of the vertical separatrix £2 has no constant part : m(Vr) = 
for all Vr- Furthermore, for the homothetic perturbation ()13p we have M(to) = for all to, 
so that the first-order Melnikov analysis is not conclusive there. In contrast, other kinds of 
perturbations, like periodical rotation around C or translation, lead to Melnikov functions of 
the form M(to) = A(uj, Vr) sin(o;io + <p) which enable to conclude that £2 is also splitted : 
this leads to extremely efficient mixing properties which could be of interest, for example, in 
microfluidic devices. Indeed, if the three separatrices break, particles released in the right-hand- 
side half-plane (x > 0) could penetrate into the cell and get mixed with particles initially located 
in the left-hand-side half-plane {x < 0). 

As noted before, the elliptic cells, like the heteroclinic bifurcation, can occur in a wide 
variety of flows. We have checked that if the basic flow is itself an elliptic dipole (like Hadamard 
cells within a bubble) with center C, then particle trajectories can take the form of a smaller 
elliptic dipole centred at C. When this basic flow is submitted to a time periodic perturbation, 
separatrices can be broken and particles can have complex trajectories within the dipole : some 
remain trapped in the Stommel cell, and some exit the cell and drop to the bubble interface. 
This model can be applied to dust transport within a drop or a bubble, a topic of interest 
in process engineering, and could help calculate bubble contamination rates. Also, for inner 
Hadamard cells or Hills' vortex one can check that 1Z = 0, so that the Taylor expansion of the 
Hamiltonian H is exact there, and the undergoing analytical results are particularly accurate. 

Finally, for non- heavy particles, buoyancy and pressure gradient and added mass can be 
readily taken into account in the particle motion equation. This leads to interesting particle 
behaviours which are currently under study. Also, Basset's drag correction, manifesting the 
effect of the unsteadiness of the disturbance flow due to the particle, could influence the particle 
dynamics and the occurrence of separatrix splitting. One can check that this correction induces 
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an (9(e 3//2 ) integral term into the asymptotic motion equation ([2]) (see also Druzhinin & Ostro- 
vsky [14J). This time integral makes the analysis much more complicated, and little results are 
available there. Note however that first-order fluid inertia effects, responsible for lift or drag 
corrections, should also be taken into account, since the Stokes number of the inclusion (mea- 
suring the unsteadiness of the particle-induced flow) is often of the order of its Taylor number 
(which characterizes first-order fluid inertia effects). Unfortunately, no general particle motion 
equation is available in this case, and the dynamics of such inclusions remains an open issue. 
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A Calculation of the Melnikov function 

The calculation of the complete integral ()10p requires to solve the differential equation 

X = V°(X (t)) + V T , (19) 
with Xq{— oo) = A and Xq{+oo) = B. Also, Xa(t) belongs to the curved separatrix, so 

a(x (t) - x c ) 2 + 7(1/0(0 - Vcf + 6 = 0, 
where Xq = (xo,yo). Therefore, the vertical coordinate of Xq satisfies : 

OH / \ 

yo = —fa= 2a(xo(t) - x c ? = -2(7(2/oW - Vcf + 6) , 

and the solution of this elementary equation, with yo(±oo) = y c + b and yo(0) = y c , is yo{t) = 
y c — b tanh(27&t) where b = \J — 6py. We readily obtain the horizontal coordinate of Xo(t) : 
xo(t) = x c + 0/ cosh(27&t), where a = yJ—5/a. The constant part of the Melnikov function is 
therefore 

f+OQ 

m= (x yo - x y ) dt 



j —00 

_4a6 3 7 2 / = -2vra6 3 7 2 . 

J-00 cosh s 



By writing a, b and 7 in terms of V c , Vt and ifP we are led to expression (|12p . 

The Melnikov integral corresponding to the homothetic perturbation (|13p . by using (|19p . reads 



M(t ) = / X AVf(X (t))smuj(t + to)dt-krn. 
Using once again (fl"9|) : 

/+00 ^ 
X A V T sinw(t + t )dt-km 
-00 



-00 

"+0O 



/+OO 
x (t)e iut dt-km, 
-00 



where Im denotes the imaginary part and Vt = —Vre y . By inserting the analytical expression 
of xo(t) into this Fourier integral we are led to : 

x (t)e lult dt = — — — , 

, 2i 07 cosh(7ra;/4o7) 

and this leads to M(to) = — A(oj, Vt) cosuto — km(Vr), with A(u, Vt) given by (fTSp . 
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B Lower bound for Vp leading to heteroclinic bifurcation 



In this appendix we make use of the Melnikov function M(to) to derive an analytical expression 
of the range of terminal velocities Vr for which the heteroclinic bifurcation does not occur. By 
using (|14p we obtain a necessary condition for the existence of both the elliptic structure and 
simple zeros in M(io) : 

V c - £>(V°) < V T < V c , (20) 

where D(ijP) is : 



s 2 a„2y _ 4 r)V3 



and D\ = Di(tp°) is given by : 



and 



£>i(^°) = 

27V;V/2 - 9V c s 4 + s e + - ^8lV c 4 s 4 - 12V c 3 s 6 (22) 

(23) 



~ D.663. Tn obtain th 

q cosh q 



and k = max — - — ~ 0.663. To obtain this result we rewrite (1141) as 



2^/6 V T 



A(u,V T ) = ^ZL^TZV^ F( = ) (24) 



with F(q) = q/ cosh One can easily check that F is a bounded function with k = max — - — 

q coshq 1 

Hence, whatever oj the function F in (|24p is at most equal to k, so that condition 

2\/6y T 



- Vt k < k \ m(Vr) 



implies that M(to) does not have any zero. By taking the square of this inequality, which 
only involves positive quantities, we are led to a condition of the form Ps(Vr) > where 
Ps(X) = {V c — X) s — s 2 X 2 is a 3rd order polynomial, and s is given by ()23p . One can easily 
check that P3 is a decaying function of X, for all X > 0, and that -Ps(O) = > and 
■f^(Vc) = — K; 2 < : P3 has a real root somewhere between and V c . If Vr is smaller than 
this root, then P^{Vt) > and M(to) does not have any zero. An analytical expression of the 
root can be found by solving P3 and we are led to Vr < V c — D(ip°) with D(ip°) given by Eqs. 
(|2ip - (|22p . This is a sufficient condition for non-chaotic sedimentation. The opposite of this 
condition provides a necessary condition for the heteroclinic tangle to occur. 

In the three dimensional axisymmetric case the same method can be applied, starting from 
the Melnikov function given by (|17p and (|18|) . We are led to the same polynomial P3, so that 
formulas ([2T]) - ([22]) can be applied, but with a different value for s : 

(25) 



q 2 

and k = max — - — ~ 1.1046. Finally, note that all these formulas are valid if 7^ = 0. 
<2>o sinhg 
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C Effect of the remainder of the streamfunction 



The particle dynamics @ can be written as : 

X p = V T + Vf(X p ) 
+ e(v}{X p ,i) - k^Vf.{Vf + V T )) + 0{e 2 ) (26) 



with 

dH 



v T + vf=[ _H 

\ dx / 



dx 

The Taylor expansion ([8]) of the particle streamfunction is of the form H(x, y) = H%(x — x c ) + 
lZ(x — x c ), where H% is a third-order polynomial function of x — x c and 1Z is a g-th order function 
of x — x c with q > 4. The inertialess particle velocity is therefore affected by the remainder : 



/ dih \ / an \ 

V dx I \ dx ' 



V 2 w i-i 

where V% is of order 2 and Wq— 1 is of order q — 1 in x — x c . The particle motion equation (|26p 

therefore reads : 

i P = V 2 + e{vj{X p ,t) - kVV 2 .V 2 ) + [W q -i 

(i) 

- ke(W 2 .W q -i + W g _i.V 2 + Wg-i.Wg-i)! (27) 

S v ' J 

(ii) 

where all the terms within the brackets are induced by 1Z and have to be negligible compared 
to the other terms. When x is in the vicinity of the ellipse we have | x — x c |= 0(max(a, b)) = 
0(| 5 | 1/2 ), since both ifi% xx (C) and i>% vy {C) in Eq. © are 0(1). Therefore, in the vicinity of 
the ellipse we have : 

W q _ x =o(| <5 \^' 2 ) < V 2 = o( I 5 I ), 
W 2 .WVi = of I <5 1 9/2 ) < W 2 .F 2 = o(\5 | 3/2 



as soon as | S |<C 1. Therefore the particle inertia term eW 2 .V 2 dominates all the terms in 
(ii). Most importantly, the particle inertia term eW 2 .V 2 also dominates the spurious term (i) 
provided ke \ 5 | 3 ^ 2 dominates | 5 J^ -1 )/ 2 , i.e. 

ke^\5\^' 2 . 

This condition is satisfied if q > 4 and if 5 is sufficiently small. (This conclusion could have also 
been obtained by normalizing (x,y) by (x/ \ 5 \ l l 2 ,y/ \ S I 1 / 2 ) in the particle motion equation.) 
In the flow ip° = sin x sin y used here we have q = 5, and this allows us to neglect the effect of the 
bracketted terms of Eq. (|27|) provided 5 is small enough. The numerical calculations presented 
here show that 5 does not need to be very small. 
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